*=============================================
* Step 4: Regressions 
// Long difference: 2008 vs. 2017
// Family FE
// IV estimates
// Note: Age is colinnear with FE --> causes error in IV estimation --> need to drop
*=============================================

use    "${dataout}MainDataset" , clear

***********************************
***** Choose years to keep!  ******
keep if year==2008 | year==2017
***********************************

keep if couple==1

/****************************************************************************/
/*******  Adjust road toll in 2017 to reflect 2008 grunnkrets ***************/
/****************************************************************************/
 
/*=== STEP 1: Get 2008-grunnkrets data =======*/
preserve
	keep if year==2008
	keep familienr indid1 indid2 grkrets grk_bed1 grk_bed2
	rename grkrets grkrets_2008 
	rename grk_bed1 grk_bed1_2008
	rename grk_bed2 grk_bed2_2008
	sort familienr indid1 indid2
	tempfile helpfile
	save `helpfile'
restore

sort familienr indid1 indid2
merge m:1 familienr indid1 indid2 using  `helpfile' , keep(match master)
drop _merge

/*=== STEP 2: Destring grunnkrets 2008 =======*/
destring grkrets_2008, gen(grkrets_2008_num)
destring grk_bed1_2008, gen(grk_bed1_2008_num)
destring grk_bed2_2008, gen(grk_bed2_2008_num)

/* === STEP 3: Adding info on 2017 road toll to 2008 grk-grk-pairs ======== */
preserve
	use "${dataout}distances", clear
	des 
	rename takst_r2017 toll_17
	keep grk1 grk2 toll_17
	tempfile grk1grk2
	save `grk1grk2'
restore

// For individual 1 and 2 in household: 
forvalues i=1(1)2 {	
	capt drop grk1
	gen grk1=grkrets_2008_num /*Use grunnkrets HOME from 2008*/

	/* to work  */
	capt drop grk2
	gen grk2=grk_bed`i'_2008_num   /*Use grunnkrets WORK from 2008*/
	merge m:1 grk1 grk2  using `grk1grk2', keep(match master)
	drop _merge
	rename toll_17 TO_toll_17_`i'

	/* changing the order of grks */
	rename grk2 temp
	rename grk1 grk2
	rename temp grk1

	/* from work */
	merge m:1 grk1 grk2  using `grk1grk2', keep(match master)
	drop _merge
	rename toll_17 FROM_toll_17_`i'

	drop grk1 grk2

	/* Creating averages of to/from work, where available */
	foreach var in toll_17_`i' {
		gen `var' = (TO_`var' + FROM_`var') / 2
		replace `var' = TO_`var' if `var' == .
		replace `var' = FROM_`var' if `var' == .
		drop TO_`var' FROM_`var'
	}
}

/* === STEP 4: Creating averages for both individuals ======== */
egen toll_17_fam_mean=rowmean(toll_17_1 toll_17_2) 

/* === STEP 5: Make adjusted road toll variable  ======== */
gen toll_fam_mean_adj_KPI=.
replace toll_fam_mean_adj_KPI=toll_fam_mean if year==2008
replace toll_fam_mean_adj_KPI=toll_17_fam_mean if year==2017  /*KPI-adjustment is using 2017 as base year - hence do not need to inflation adjust road toll*/

* ==============================================================================
// Run regressions
* ==============================================================================

// Reduced form and IV estimation
local treatvar1 toll_fam_mean
local IV1 toll_fam_mean_adj_KPI

** Select variables to loop over: 
local outcomelist ICE_fam_count   /*BEV_fam_yes ICE_fam_count car_fam_count car_fam_yes  */               
local treatvar2list ptl_fam_km_mean     

** Regressions - loop over: outcome, treatment, year
foreach outcome in `outcomelist'{
	foreach treatvar2 in `treatvar2list'{		
		// (a): IV estimation ==================================================	
		mata: mata mlib index   /*Need to include this to make ivreghdfe work*/

		ivreghdfe `outcome'  `treatvar2'  /*$age*/   $distance  $time  $publictime  $publicquality  (`treatvar1' = `IV1') , ///
		absorb($FE   $household  $income  $employment  $education  i.year  i.familienr_num ) cluster($clustervar) ///
		first ffirst savefirst savefprefix(first_) rf saverf saverfprefix(reduced_) 
		/**/
		summarize `outcome' if e(sample)==1
		estadd scalar MeanDep = r(mean), replace
		local MeanDep = r(mean)  
		summarize `treatvar1' if e(sample)==1
		estadd scalar Treat1 = r(mean), replace
		local Treat1 = r(mean)  
		summarize `IV1' if e(sample)==1
		estadd scalar Treat1IV = r(mean), replace
		local Treat1IV = r(mean) 
		summarize `treatvar2' if e(sample)==1
		estadd scalar Treat2 = r(mean), replace
		local Treat2 = r(mean) 
		/**/
		summarize `outcome' if e(sample)==1 & `treatvar1'==0 & `treatvar2'==0  
		estadd scalar MeanDepCon = r(mean), replace  
		local MeanDepCon = r(mean)  
		/**/
		local Proportional1 = _b[`treatvar1']/`MeanDepCon'  
		estadd scalar Proportional1 =  `Proportional1' 
		local Elasticity1 = (_b[`treatvar1']/`MeanDep') / (1/`Treat1')   
		estadd scalar Elasticity1 =  `Elasticity1'
		/**/
		local Proportional2 = _b[`treatvar2']/`MeanDepCon'  
		estadd scalar Proportional2 =  `Proportional2'
		local Elasticity2 = (_b[`treatvar2']/`MeanDep') / (1/`Treat2')  
		estadd scalar Elasticity2 =  `Elasticity2'
		/**/
		estadd local year="2008,2017" , replace
		estadd local grkFE="\checkmark", replace
		estadd local grkbFE="\checkmark", replace
		estadd local HHcon="\checkmark", replace
		estadd local HHWork="\checkmark", replace
		estadd local HHFE="\checkmark", replace
		estadd local PublicTime="\checkmark", replace
		eststo tableC5reg3 

		estimates save "${ster}TableC5_SecondStage" , replace

		** Mark main sample
		capt drop MainSample
		gen MainSample=0
		replace MainSample=1 if e(sample)==1

		// First stage =========================================================
		estimates restore first_`treatvar1'

		estadd local year="2008,2017" , replace
		estadd local grkFE="\checkmark", replace
		estadd local grkbFE="\checkmark", replace
		estadd local HHcon="\checkmark", replace
		estadd local HHWork="\checkmark", replace
		estadd local HHFE="\checkmark", replace
		estadd local PublicTime="\checkmark", replace

		estimates save "${ster}TableC5_FirstStage" , replace

		// Reduced form ========================================================
		estimates restore reduced_`outcome'

		summarize `outcome' if e(sample)==1
		estadd scalar MeanDep = r(mean), replace
		local MeanDep = r(mean)  
		summarize `treatvar1' if e(sample)==1
		estadd scalar Treat1 = r(mean), replace
		local Treat1 = r(mean)  
		summarize `IV1' if e(sample)==1
		estadd scalar Treat1IV = r(mean), replace
		local Treat1IV = r(mean)
		summarize `treatvar2' if e(sample)==1
		estadd scalar Treat2 = r(mean), replace
		local Treat2 = r(mean)  
		/**/
		summarize `outcome' if e(sample)==1 & `IV1'==0 & `treatvar2'==0  
		estadd scalar MeanDepCon = r(mean), replace  
		local MeanDepCon = r(mean) 
		/**/
		local Proportional1 = _b[`IV1']/`MeanDepCon'  
		estadd scalar Proportional1 =  `Proportional1' 
		local Elasticity1 = (_b[`IV1']/`MeanDep') / (1/`Treat1IV')   
		estadd scalar Elasticity1 =  `Elasticity1' 
		/**/
		local Proportional2 = _b[`treatvar2']/`MeanDepCon' 
		estadd scalar Proportional2 =  `Proportional2'
		local Elasticity2 = (_b[`treatvar2']/`MeanDep') / (1/`Treat2')  
		estadd scalar Elasticity2 =  `Elasticity2'
		/**/
		estadd local year="2008,2017" , replace
		estadd local grkFE="\checkmark", replace
		estadd local grkbFE="\checkmark", replace
		estadd local HHcon="\checkmark", replace
		estadd local HHWork="\checkmark", replace
		estadd local HHFE="\checkmark", replace
		estadd local PublicTime="\checkmark", replace

		eststo reduced_`outcome'
		estimates save "${ster}TableC5_ReducedForm" , replace

		// (b): "Normal"  OLS estimation with toll_fam_mean_KPI ================
		reghdfe `outcome' `treatvar1'  `treatvar2'  /*$age*/    $distance  $time  $publictime  $publicquality  if MainSample==1 , ///
		absorb($FE   $household  $income  $employment  $education  i.year  i.familienr_num ) vce(cluster $clustervar) 
		/**/
		summarize `outcome' if e(sample)==1
		estadd scalar MeanDep = r(mean), replace
		local MeanDep = r(mean) 
		summarize `treatvar1' if e(sample)==1
		estadd scalar Treat1 = r(mean), replace
		local Treat1 = r(mean) 
		summarize `IV1' if e(sample)==1
		estadd scalar Treat1IV = r(mean), replace
		local Treat1IV = r(mean) 
		summarize `treatvar2' if e(sample)==1
		estadd scalar Treat2 = r(mean), replace
		local Treat2 = r(mean)
		/**/
		summarize `outcome' if e(sample)==1 & `treatvar1'==0 & `treatvar2'==0 
		estadd scalar MeanDepCon = r(mean), replace   
		local MeanDepCon = r(mean)  
		/**/
		local Proportional1 = _b[`treatvar1']/`MeanDepCon'  
		estadd scalar Proportional1 =  `Proportional1'
		local Elasticity1 = (_b[`treatvar1']/`MeanDep') / (1/`Treat1')  
		estadd scalar Elasticity1 =  `Elasticity1' 
		/**/
		local Proportional2 = _b[`treatvar2']/`MeanDepCon'  
		estadd scalar Proportional2 =  `Proportional2' 
		local Elasticity2 = (_b[`treatvar2']/`MeanDep') / (1/`Treat2')   
		estadd scalar Elasticity2 =  `Elasticity2'
		/**/

		estadd local year="2008,2017" , replace
		estadd local grkFE="\checkmark", replace
		estadd local grkbFE="\checkmark", replace
		estadd local HHcon="\checkmark", replace
		estadd local HHWork="\checkmark", replace
		estadd local HHFE="\checkmark", replace
		estadd local PublicTime="\checkmark", replace
		eststo tableC5reg2
		estimates save "${ster}TableC5_FamilyFE_OLS" , replace

		// (c): "Normal"  OLS estimation with toll_fam_mean_KPI and without family FE 
		reghdfe `outcome' `treatvar1'  `treatvar2'  /*$age*/    $distance  $time  $publictime  $publicquality  if MainSample==1 , ///
		absorb($FE   $household  $income  $employment  $education  i.year   ) vce(cluster $clustervar) 
		/**/
		summarize `outcome' if e(sample)==1
		estadd scalar MeanDep = r(mean), replace
		local MeanDep = r(mean) 

		summarize `treatvar1' if e(sample)==1
		estadd scalar Treat1 = r(mean), replace
		local Treat1 = r(mean)  

		summarize `IV1' if e(sample)==1
		estadd scalar Treat1IV = r(mean), replace
		local Treat1IV = r(mean)  

		summarize `treatvar2' if e(sample)==1
		estadd scalar Treat2 = r(mean), replace
		local Treat2 = r(mean)  
		/**/
		summarize `outcome' if e(sample)==1 & `treatvar1'==0 & `treatvar2'==0  
		estadd scalar MeanDepCon = r(mean), replace  
		local MeanDepCon = r(mean) 
		/**/

		local Proportional1 = _b[`treatvar1']/`MeanDepCon'   
		estadd scalar Proportional1 =  `Proportional1' 
		local Elasticity1 = (_b[`treatvar1']/`MeanDep') / (1/`Treat1')   
		estadd scalar Elasticity1 =  `Elasticity1' 
		/**/
		local Proportional2 = _b[`treatvar2']/`MeanDepCon'   
		estadd scalar Proportional2 =  `Proportional2' 
		local Elasticity2 = (_b[`treatvar2']/`MeanDep') / (1/`Treat2')  
		estadd scalar Elasticity2 =  `Elasticity2' 
		/**/
		estadd local year="2008,2017" , replace
		estadd local grkFE="\checkmark", replace
		estadd local grkbFE="\checkmark", replace
		estadd local HHcon="\checkmark", replace
		estadd local HHWork="\checkmark", replace
		estadd local HHFE="", replace
		estadd local PublicTime="\checkmark", replace
		eststo tableC5reg1 
		estimates save "${ster}TableC5_GrunnkretsFE_OLS" , replace
	}
}

*=================================================
*=================================================
*===== Make table from estimates saved to disk
*=================================================
*=================================================

** Select variables to loop over: 
local outcomelist   ICE_fam_count   /*BEV_fam_yes ICE_fam_count*/
local treatvar1list toll_fam_mean     
local treatvar2list ptl_fam_km_mean    
local IV1list toll_fam_mean_adj_KPI     

foreach outcome in `outcomelist'{
	foreach treatvar1 in `treatvar1list'{
		foreach treatvar2 in `treatvar2list'{
			foreach IV1 in `IV1list'{				
				eststo drop *
				* Load estimates
				estimates use "${ster}TableC5_GrunnkretsFE_OLS"
				eststo
				estimates use "${ster}TableC5_FamilyFE_OLS" 
				eststo
				estimates use "${ster}TableC5_SecondStage" 
				eststo
				estimates use "${ster}TableC5_FirstStage" 
				estadd local year="2008,2017" , replace
				estadd local grkFE="\checkmark", replace
				estadd local grkbFE="\checkmark", replace
				estadd local HHcon="\checkmark", replace
				estadd local HHWork="\checkmark", replace
				estadd local HHFE="\checkmark", replace
				estadd local PublicTime="\checkmark", replace
				eststo
				estimates use "${ster}TableC5_ReducedForm" 
				estadd local MeanDep="" , replace
				estadd local Treat1="" , replace
				estadd local Treat1IV="" , replace
				estadd local Proportional1="" , replace
				estadd local Elasticity1="" , replace
				eststo

				// Latex 
				esttab est*  ///
				using  "${tables}TableC5.tex" ///
				, keep(`treatvar1'  `IV1') ///
				title("Robustness: Effects of road toll on BEV ownership, 2008 to 2017 ") ///
				booktabs nonotes  /*nomtitle  label*/  coeflabels(`treatvar1' "Road toll (NOK)"  `IV1' "Road toll 2008 commute (NOK)") ///
				mtitles("OLS" "OLS" "\shortstack{Second \\stage}"  "\shortstack{First \\stage}"  "\shortstack{Reduced \\form}") ///
				stats( N year grkFE grkbFE HHFE  HHWork HHcon  PublicTime MeanDep Treat1  Treat1IV widstat Proportional1 Elasticity1 , label("N" "Year"  "Neighborhood residence FE" "Neighborhood work FE"  "Household FE"  "Ind-work controls" "HH and ind controls" "Public transit (minutes)" "Mean outcome" "Mean road toll (NOK)"  "Mean road toll 2008 commute (NOK)" "F statistic (excl. instrument)" "Road toll, prop. effect" "Road toll elasticity" ) fmt(%9.0fc a3) ) ///
				b(5) se(5) star(* 0.10 ** 0.05 *** 0.01)    replace

				// Latex  - without table opening
				esttab est*  ///
				using  "${tables}TableC5Fragment.tex" ///
				, keep(`treatvar1'  `IV1') ///
				title("Robustness: Effects of road toll on BEV ownership, 2008 to 2017 ") ///
				booktabs nonotes  /*nomtitle  label*/  coeflabels(`treatvar1'  "Road toll (NOK)"   `IV1' "Road toll 2008 commute (NOK)") ///
				mtitles("OLS" "OLS" "\shortstack{Second \\stage}"  "\shortstack{First \\stage}"  "\shortstack{Reduced \\form}") ///
				stats( N year grkFE grkbFE HHFE  HHWork HHcon  PublicTime MeanDep Treat1  Treat1IV widstat Proportional1 Elasticity1 , label("N" "Year"  "Neighborhood residence FE" "Neighborhood work FE"  "Household FE"  "Ind-work controls" "HH and ind controls" "Public transit (minutes)" "Mean outcome" "Mean road toll (NOK)"  "Mean road toll 2008 commute (NOK)" "F statistic (excl. instrument)" "Road toll, prop. effect" "Road toll elasticity" ) fmt(%9.0fc a3) ) ///
				b(5) se(5) star(* 0.10 ** 0.05 *** 0.01) fragment   replace

				// To screen
				esttab est*  ///
				, keep(`treatvar1'  `IV1') ///
				title("Robustness: Effects of road toll on BEV ownership, 2008 to 2017 ") ///
				booktabs nonotes  /*nomtitle  label*/  coeflabels(`treatvar1'  "Road toll (NOK)"   `IV1' "Road toll 2008 commute (NOK)") ///
				mtitles("OLS" "OLS" "\shortstack{Second \\stage}"  "\shortstack{First \\stage}"  "\shortstack{Reduced \\form}") ///
				stats( N year grkFE grkbFE HHFE  HHWork HHcon  PublicTime MeanDep Treat1  Treat1IV widstat Proportional1 Elasticity1 , label("N" "Year"  "Neighborhood residence FE" "Neighborhood work FE"  "Household FE"  "Ind-work controls" "HH and ind controls" "Public transit (minutes)" "Mean outcome" "Mean road toll (NOK)"  "Mean road toll 2008 commute (NOK)" "F statistic (excl. instrument)" "Road toll, prop. effect" "Road toll elasticity" ) fmt(%9.0fc a3) ) ///
				b(5) se(5) star(* 0.10 ** 0.05 *** 0.01)

			}
		}
	}
}